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Abstract 



Starting from considerations about meaning and subsequent use of asymmetric 
uncertainty intervals of experimental results, we review the issue of uncertainty 
propagation. We show that, using a probabilistic approach (the so-called Bayesian 
approach), all sources of uncertainty can be included in a logically consistent way. 
Practical formulae for the first moments of the probability distribution are derived 
up to second-order approximations. 
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1 Introduction 

The combination in quadrature of uncertainties due to systematic effects has become quite 
standard practice in physics. It is also common practice to add these uncertainties in 
quadrature to those from random effects. Usually the two kinds of uncertainties are given 
separately, and the systematic-effect uncertainties are listed individually (at least for the 
most relevant ones) in order to show the potential of further measurements made with the 
same apparatus. This combination rule has arisen as a kind of pragmatic procedure |lj , in 
analogy to the combination of standard deviations in probability theory, although cannot 
justifiably be termed within 'conventional' (i.e. non-Bayesian) statistics. The same is true 
for the use of the covariance matrix to handle correlated uncertainties. 

There is less agreement when the uncertainties due to systematic effects are asym- 
metric and/or they produce asymmetric shifts in the final quantity of interest, due to 
nonlinear propagation of uncertainty.^] As a numerical example of the latter case, take 
a quantity Y depending on three 'influence parameters' X\, X 2 and X 3 , which could be 
calibration constants, environment quantities or theoretical parameters. Suppose that, 
for the reference values of the X's, the analysis procedure gives (in arbitrary units) 
Y = 1. 000 ± 0.050, where the uncertainty associated with the result is that due to random 
effects. Consider now that by 'varying reasonably the parameters Xf (the expression is 
intentionally left vague for the moment) the following deviations from the central values 
occur: AY l± = tomo, AF 2 ± = tai47> and AY 3± = ±ai56- An often-used practice is to 
combine in quadrature separately positive and negative deviations, obtaining the follow- 
ing result: Y = 1.00 ± 0.05 (stat.) to'^t ( S Y S ^-)- Now we are faced with the problem that 
the result of this ad hoc procedure has no theoretical justification. Hence the uncertainty 
content of the statement (i.e. its probabilistic meaning) is unclear and, as a consequence, 
it is not obvious how to make use of this information in further analyses, even in the 
simple case in which the data points are uncorrected. As a matter of fact, most people 
remove the asymmetry in further analysis of the results, using something equivalent to 
the standard deviations to be used in x 2 fits. This 'standard deviation' is evaluated either 
by taking the largest value between A + and A_, or by averaging the two values (some 
use the arithmetic, others the geometric average). The result is that in both procedures 
the uncertainty is symmetrized and the result is considered as if it were described, for 
all practical purposes, by a Gaussian model around the published best estimate^] Our 
main worry is not that the combined uncertainties will be incorrect (we anticipate that 

^ As is well known, asymmetric uncertainties arise also from random effects alone, for example in \ 2 fits 
if the x 2 around the minimum is not symmetric. In this paper we focus only on asymmetries deriving 
from systematic effects or non-linearity. 

2 ) A different problem, although also related to systematic effects, is the treatment of overall uncertainty 
common to all data, or to large subsets of data, in fitting. In fact, it is now understood that the 
covariance matrix techniques might lead to undesirable effects on the results We intend to 
review this problem in a separate paper. 

3 ) A more complicated 'prescription' is described by the PDG ||, which we report here for the con- 
venience of the reader: "When experimenters quote asymmetric errors (Sx) + and (6x)~ for a mea- 
surement x, the error that we use for that measurement in making an average or a fit with other 
measurements is a continuous function of these three quantities. When the resultant average or fit x 
is less than x— (5x)~ , we use (Sx)~ ; when it is greater than x+ (Sx) + , we use (5x) + . In between, the 
error we use is a linear function of x. Since the errors we use are functions of the result, we iterate to 
get the final result. Asymmetric output errors are determined from the input errors assuming a linear 
relation between the input and the output quantities. " This rule does not seem to be applied by others 
then the PDG. As examples of other ad hoc procedures, see Refs. [pHq]. 
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the arithmetic average of A + and A_ gives indeed the correct uncertainty in most cases 
of practical interest), but rather that the result itself can be biased with respect to what 
one could get using consistently the best knowledge concerning the input quantities, as 
will be shown in this paper. 

The purpose of this paper is to review the issue of uncertainty propagation, starting 
from general considerations and deducing approximate formulae for practical applications. 
The issue will be analysed in the framework of the so-called Bayesian inference, though 
Bayes' theorem will never appear in this paper. This just means that, following physics 
intuition we consider it natural to talk about probability of true values. As a conse- 
quence, for this kind of application probability can only have the meaning of degree of 
belief.0 We will show that the rule of combining in quadrature symmetric uncertainties is 
a natural consequence of the probabilistic approach we follow, assuming that the uncer- 
tainties are conceptually properly defined (although the overall result will not depend on 
their precise values). Formulae to take into account correlations and nonlinearity effects 
will also be provided. 

The paper is structured as follows. First, in Section ^| we illustrate briefly what we 
mean by probabilistic treatment of measurement uncertainties and why this is meaning- 
ful only within the Bayesian approach. In this approach, the most general and logically 
consistent way to include all kinds of uncertainties is just a straightforward application 
of probability calculus. However, the application of the general propagation formula (the 
derivation of which is given in Appendix A) requires the evaluation of integrals which 
can rapidly become complicated in real-life problems. Therefore, approximate formulae 
are derived based on linear and quadratic expansion of the output quantity on the input 
quantities. This is done in Sections |3| and [5], respectively. Several numerical examples are 
given, the main ones being discussed in Section |6|. Section |7| shows that our approximate 
formulae can also handle the case in which ±1 a variations on an input quantity produce 
a shift in the same direction on the output quantity. Correlations are also considered. 
We describe in the text only the linear case (Section |3|) and refer to Appendix B for the 
more complicated formulae which take into account second-order effects. The important 
issue of how to model uncertainties due to systematic effects is discussed in Section f|. For 
practical purposes the probability density function (p.d.f.) of the output quantity can be 
considered approximately Gaussian also in the non-linear cases and also if correlations are 
present, thanks to the combinatorial effects analogous to those which make the central 
limit theorem work. In case of important non-Gaussian contributions in the input quan- 
tities, or strong non-linearity effects in the propagation, a detailed evaluation of the p.d.f. 
is needed, usually done using Monte Carlo methods. However, deviations from normality 
can be checked from skewness and kurtosis, and we give approximate formulae for these 
quantities in Appendix B. Finally, we show in Appendix C how it is possible, in principle, 
to get a rough evaluation of the final p.d.f., if this does not differ much from a Gaussian; 
'in principle' because we understand that the method described in Appendix C is perhaps 
more an academic exercise than a real help to practitioners, who most likely will find it 
more convenient to solve the problem by Monte Carlo integration. Concluding remarks 
are given in Section 



4 ) For a physicist's introduction to subjective probability and Bayesian inference see Ref. ||, or Refs. jio] 
and pin for short accounts. 
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2 Probabilistic treatment of measurement uncertainties 



In the Bayesian approach, probability is associated with uncertainty, whether we are 
interested in the yet to be observed outcome of a measurement, or in the numerical value 
of the physics quantity. As a result of the experiment, there will be a p.d.f. /(/x | data, I) 
associated with the numerical value of the true value (generically called /i), conditioned 
by the observed data and by the status of information (/) concerning measurement and 
measurand (the conditions data and / will usually be considered implicit). Although we 
maintain that the proper way of learning from data is to make use of Bayes' theorem, 
it is easy to show that in most routineQ (or, at least, non critical) measurements the 
usual methods of analysis can be considered as approximations of Bayesian inference 
(see also Section 2.9 of Ref. ||). When these kinds of conditions hold, also the Gaussian 
approximation is usually rather good. Therefore, hereafter we will consider a result of the 
kind /i = ft ± a r , where a r is only due to random effects, equivalent to a Gaussian p.d.f. 
of /i, as usually perceived by physicists |],f|], i.e. 



(Mr -A) 2 



/i = ft ± cr r f(n r I data, I) = — — e 2a ? , (1) 

V 2 7T o r 

where the symbol /i r is to remind us that only uncertainties due to random effects are 
considered in Eq. (|l|). On the other hand, if one sticks strictly to frequentistic ideas, 
one gets 'results' which have neither a meaning of certain statements about true values, 
neither the meaning of probabilistic statements. As a consequence, it is not clear what 
they mean, neither how they should be combined or propagated in a logically consistent 



way (see more extended discussion in Refs. |13| and |L4 



In the probabilistic framework in which we are moving, the uncertainties due to 
systematic effects can be easily and consistently included (at least conceptually, although 
the numerical implementation can present some technical problems). Indeed, there are 
several ways to proceed, all leading to the same result, though each way can be more or 
less intuitive or suitable for a particular application (see Section 2.10.3 of Ref. ||). The 
closest one to the spirit of probabilistic inference consists in writing explicitly I to depend 
on other physics quantities, which could be calibration constants, influence parameters 
(temperature, pressure, etc), theoretical quantities, and so on, plus other pieces of general 
knowledge not easy to model (I ) and which lead the researchers to behave in a given 
way and to make reasonable assumptions in the many steps of the experimental work. 
Let us indicate by I this general knowledge. The physical quantities on which the result 
can depend will be called influence quantities (or parameters) and will be indicated by hi. 
The entire set of influence quantities will be indicated by h = {hi, hi, ■ ■ ■ , h n }. In general, 
the result ([I]), which takes into account only random effects, is obtained using the set of 
best estimates^] of the parameters (h Q ), i.e. 

/(/i r | data, /) = f(fi | data, h D , I Q ). (2) 

5 ) 'Routine' in the sense of Ref. |L2]|, which applies also for most measurements in frontier research. 

6 ) Note that sometimes the result is meant, instead, for a nominal set of parameters, which are not 
necessarily the ones in which the experimentalists believe mostly. In principle, it is possible to derive 
the approximate formulae considering nominal values, as has been done by one of the authors in 
Ref. |^| . However, when second-order effects are taken into account the formalism becomes more 
complicated and, therefore, we prefer to start from expected values and standard deviations on the 
influence quantities. The reader should be aware that the use of nominal values different from expected 
values produce shifts in the results which need to be corrected similar to those produced by 
nonlinearity effects which will be discussed later in this paper. 
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The most general inference on fi will depend, instead, on all possible values of h, and 
the resulting p.d.f. will be f(fi | data, I Q ). Probability theory teaches us how to get rid 
of the uncertainty about the exact value of the influence parameters. Describing the 
uncertainty about the influence parameters with the joint p.d.f. f(h | I Q ), we obtain that 
the probabilistic result which takes into account systematic uncertainties is given by 

/(/i | data, I Q ) = J fi.^ I data, h, J Q ) f(h | J D ) dh . (3) 

(We use the symbol /(•) for all p.d.f.'s, and implicitly consider the integrals done over 
the range of definition of the variables.) If the influence parameters are perfectly known, 
i.e. f(h 1 1 ) = Ui5(hi — h 0i ), we reobtain Eq. (H), and hence Eq. (|l|), i.e. the uncertainty 
is that due to the random effects alone. Hereafter we shall consider implicit the general 
condition I Q . 

As an example, let us consider the result of a single measurement yielding the 
observed value X = x, and in which the most relevant systematic effect is a not exactly 
known offset Z, the uncertainty about which is described by a Gaussian p.d.f. around zero 
and standard deviation a z . We have (see Ref. [|1(J for further details): 

1 _(v-x-z) 2 

f(fji\x,z) = e , (4) 

/ + °° l _( M -*-z) 2 X _ z 2 

i= e e ^dz (5) 

-oa y In o r y in a z 

1 (m-z) 2 

= ■ ^^ =e n<^+«i) . (6) 

The p.d.f. which describes \i is still centred around the observed value x, but with a 
standard deviation which is the quadratic combination of o r and a z . The commonly used 
combination rule is recovered, but now as a theorem with well-defined conditions, instead 
of just a 'prescription'. 

An alternative way of including systematic effects, very convenient for deriving 
approximate formulae, consists in considering a function g which relates the true value \x 
to \i T and of the influence factors, i.e. 

/i = g(/i r ,h) . (7) 

Therefore the uncertainty about \l is obtained from the propagation of uncertainties about 
/i r and h (see Appendix A): 



f(p) = J /(A*r,h) ■ 5( / u-5'(/ir,h))d/i r dh. 



This formula also has a simple interpretation which makes it convenient for Monte Carlo 
evaluation^ the infinitesimal probability element /(//) d/^ depends on 'how many' ele- 
ments d/i r dh contribute to it, each element weighted with the p.d.f. calculated in {/i r , h}. 

7 ) For example, this is the basic reasoning behind the methods used by several authors to evaluate the 
p.d.f.'s of physical quantities, like the direct CP-violation parameter e'/e Jl5[|l6|, or the parameters 
of the quark mixing matrix |]l6| , : Beliefs on the input quantities (experimental and theoretical 
quantities) are propagated into the beliefs on e'/e, or on p and ij, respectively. The result has a clear 
probabilistic meaning and is, as we shall see, rather insensitive to the exact shape of the input p.d.f.'s. 
Instead, the so-called 'scanning' |l6| or other ad hoc procedures (see e.g. Ref. |i"8[| ) do not have 
such an intuitive interpretation and can be misleading, especially when the very conservative 'regions 
of confidence' produced by these methods are improperly called 95% C.L. regions. 
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At this point, an interesting observation is that fi r and hi have a symmetric role in 
the propagation of uncertainty, and therefore there is no real need to keep them separate 
in the formalism. Therefore, following the ISO Guide | ]L9|| , we prefer to speak, generically, 
of input quantities, and to indicate them all by Xj. We indicate the output quantity by Y . 
In many problems of interest also the output quantities might be more than one. Their 
values are evaluated using the common data or (which is conceptually equivalent) making 
use of the same instrumentation. In such a case we have to consider correlations among 
the output quantities even if the input quantities were uncorrelated. Hereafter we will 
indicate the generic functions gi(X) with the same symbol as the output variables, and 
speak about Yj = 3^(X). 



3 Linear expansion around E[X], role of central limit theorem and 
numerical implementation of the linear propagation 

Having illustrated the general solution to the problem, it is now interesting to obtain 
approximate formulae which, in many practical cases, save us from making complicated 
integrals. The case in which the dependence of Yj on Xi is approximately linear in a range 
of several standard deviations around their expected value is well known, and leads to the 
standard propagation formula of variances and covariances. What is less well known is that 
the use of these formulae is justified only if the numerical values of the physics quantities 
are associated with random numbers (or uncertain numbers), and the probability is meant 
as degree of belief |9|,|l0fl. 



The first-order expansion of (X) around the expected values of Xi gives 
Yj « Y}(E[X])+^ 

Xi, (10) 



dY^ 
dXi 



(Xi-E[Xi\) (9) 

E[x] 



dYj 

dX. 



E[X] 



where E[-] stands for expected value and the derivatives are evaluated for x = E[X] (this 
will be implicit hereafter). The second formula is very convenient to calculate the variance, 
having put in k all terms which do not contain Xi. Evaluating the expected values from 
Eq. (HO , and variances and covariances from Eq. ( |10|) , we get (we have replaced the symbol 
by '=' to indicate that there are no further approximations than linearization): 

E [ Y j] = ^(E[X]), (11) 

- 2 (^) = >.(^) "f-i m^-i^} ■ (12) 




where <7j are shorthand for erpQ) and pi m are the correlation coefficients, such that 
Pim°~io~ m = Cov(Xi,X m ). The terms within {•} vanish if the input quantities are un- 
correlated, as it often the case when relevant systematic effects are considered. However, 
sometimes this is not the CclSG, clS when several calibration constants are simultaneously 
obtained from a fit. Equations (|1^)-(|13D can be written in the more compact form of 
covariance matrix transformation. However, for the purpose of this paper, we prefer not 
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to use the matrix formalism, in order to separate the contributions due to variances and 
covariances. 

Equations ([TTD (P~5|) give only some limited information about the joint p.d.f. of Y, 
namely only 1st and 2nd moments. However the central limit theorem plays the important 
role of making the p.d.f. of each Yj practically Gaussian in most of the cases of interest 
(see e.g. examples and words of caution in Ref. 0, and discussion in Ref. [pi]] ). The joint 
p.d.f. can be considered for practical purposes a multivariate Gaussian. Anyway, in case 
of doubt, it is good practice to check the shape of each marginal p.d.f. (see Appendix C). 

In complex real-life cases the derivatives are not performed analytically. Instead, 
the effects of the input values on the output values are evaluated numerically, often by 
Monte Carlo techniques. In these cases the derivatives can be estimated numerically by 
± 1 a variations around the expected values. Calling A± the variation of Yj due to a 
variation of Xi of ±1 (Tj around E[Aj], linearity implies that 



OA, Oi (Ji 



Since in the linear approximation A +ji and A . 4 are practically equal, we call Aji either 
of them (taking the average of the two if there are small differences; the case of large 
differences, hint of non-linear effects, will be discussed below). We get, finally, the following 
practical formulae for the elements of the covariance matrix: 



a 2 (Yj) = + \ 2 E P>™ A ^ } > 

i \ Km J 

i \ l<m 



(15) 

Cov(Y}, Y k ) = V Aji A kl + I 2 V pi m A j7 A km } . (16) 



In the simple case of independent input quantities, Eqs. ([I5|)- (ll6|) reduce to 

a\Yj) = -£A% (17) 



Cov(y i ,y fe ) = J2A j{ A ki 



22 Cov^Yj, Y k ) = 2j s ijk \Aji\ \A ki \ 



where CoVj(Yj, Y k ) stands for the contribution to the covariance from the zth input quan- 
tity, and Sij k indicate the product of the signs of the absolute increments of Yj and Y k for 
a variation of Xi (|Ajj| have the meaning of standard uncertainty of Yj due to Xi alone). 

At this point, we have to remember that fi r defined in Sec. |] is considered as one 
of the input quantities, and that in the most general case there will be many \x T ., each 
associated with one and only one output quantity Yj. The resulting covariance matrix 
will be equal to the sum of the covariance matrix of the /i r . (they can be correlated as 
they could come from fitting procedures, unfolding, or other statistical techniques) and 
the covariance matrix due to the systematic effects. Let us write down, as an easy and 
practical example, the formulae for the case when we have iV values /v and the influence 



We have used the following notation: A + = Y(E[X] + a x ) - ^(E[X]) and A_ = Y(E[X])- Y(E[X] - 
<tx)- Therefore, for monotonic functions around E[X] the increments A+ and A_ have the same sign. 
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quantities are uncorrelated: 



^) = <+E A ?- ( 19 ) 

i>N 

Cov(Yj, Y k ) = Cov(n r .,n r .) + ^ s ijk \^ki\, (20) 

i>N 

where we have taken into account that the Aji associated with fi n are given by Aji = <Ji Sij, 
where 5^ is Kronecker symbol. In fact, the derivatives of Yj with respect to [i r ., evaluated 
at the point of best estimate of X, are equal to 1 if i = j, and equal to otherwise. 

4 Modelling the uncertainty due to systematic effects: ISO type-B 
uncertainties 

At this point it is important to define somewhat better how the several ingredients ap- 
pearing in the previous formulae should be evaluated. In fact, the results of the above 
formulae have a defined probabilistic meaning only if the various A's are obtained as 
variations of the output quantities for 1 a variations of the input quantities, and not, 
generically, as reasonable variations, or, prudentially, as 'conservative variations'. Now we 
are confronted with the problem that in the evaluation of uncertainties due to imperfect 
knowledge of systematic effects, the case in which the input uncertainties are evaluated 
from standard statistical procedures which provide standard deviations in an automatic 
way is rare. These latter cases would be those in which we feel comfortable. More often, 
"for estimate Xi of an input quantity Xj that has not been obtained from repeated obser- 
vations, the . . . standard uncertainty . . . is evaluated by scientific judgement based on all 
the available information on the possible variability of Xi. The pool of information may 
include: previous data; experience with or general knowledge of the behaviour and prop- 
erties of relevant materials and instruments; manufacturer's specifications; data provided 
in calibration and other certificates; uncertainties assigned to reference data taken from 
handbooks." (ISO Guide 0]). This is along the spirit that "the evaluation of uncertainty 
is neither a routine task nor a purely mathematical one; it depends on detailed knowledge 
of the nature of the measurand and of the measurement" |L9] . 



Following the recommendation of the BIPM recommendations ||20|1 , the ISO Guide 
calls this kind of uncertainty type B, in contrast to type-A uncertainties obtained, to say 
it shortly, by the dispersion of readings (see also Section 6.1.2 of Ref. |§). The evaluation 
of type-B uncertainties implies the adoption of the "viewpoint . . . that probability is a 
measure of the degree of belief that an event will occur" [[0|]. In practice it requires a 
realistic and honest modelling of the case. The most common models are discussed in 
the ISO Guide itself: For example, if one is practically sure that an input value is in a 
certain interval, and all values inside the interval appear similarly likely, the proper model 
for the uncertainty is a uniform distribution. Other times the edges of the interval seem 
still to be really extreme values for the quantity; but one tends to believe more in central 
values, and the belief decreases roughly linearly from the centre to the edges. In this 
case, a more suitable distribution is a symmetric triangular distribution. Alternatively, if 
the belief decreases towards the edges, but the maximum belief does not coincide with 
(approximately) the centre of the interval, it is preferable to use an asymmetric triangular 
distribution. Finally, if the interval seems simply to be just a highly probable one (e.g. 90%, 
95% or 99%), but also far away values are believed to be possible, one can use a Gaussian 
model with a standard deviation which is a suitable fraction of the uncertainty interval. 
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Figure 1: Typical models to assess type-B uncertainties: uniform distribution, symmetric trian- 
gular distribution, asymmetric triangular distribution, and Gaussian distribution. The expres- 
sions of the most relevant statistical parameters are reported (S stands for skewness, K, for 
kurtosis). 



Figure [I] shows the most common models to handle type-B uncertainties, together with 
their most interesting statistical parameters. 

Once the mathematical model has been chosen, the p.d.f. of the output quantity can 
be evaluated using Eq. or their first moments can be obtained using approximated 
formulae. It is important to realize that neither the choice of model, nor the value of the 
standard deviation are very critical, as discussed in Refs. and |19| (see also Section 4.6 
of Ref. |21|| ). In fact, the central limit theorem makes the result Gaussian independently 
of the initial distribution, if none of the non-Gaussian components dominates. Moreover, 
some unavoidable over- or under-estimates, compensate, if one makes the effort of assessing 
model and standard deviation in an honest way. Finally, in agreement with the ISO Guide, 
we think that it is bad practice to overestimate intentionally type-B uncertainties (see 
Section 10.2 of Ref. §). 
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As a numerical example, let us consider the standard deviations of input quantities 
believed to be, with certainty or with high probability, in the interval between —1 and 
+1. 

Uniform: cr(X) 
Symmetric triangular: &{X) 

Asymmetric triangular peaked at 1/2: &{X) 
Gaussian, 90% probability interval: &{X) 
Gaussian, 95% probability interval: cr(X) 

We see that, for practical purposes, the differences between the ex's are irrelevant. Never- 
theless, in order to avoid a bias of the overall uncertainty, one should try to model each 
component according to the best knowledge of the physics case, rather than by choosing 
systematically the model which gives the most conservative uncertainty^Note that in the 
case of asymmetric triangular distribution, the expected value of X is neither the centre 
of the interval, nor the peak of the distribution. In this case we have E[X] = 1/6 ~ 0.17. 
If one uses, incorrectly, the peak value, one introduces a bias which is ~ 80% of a stan- 
dard deviation. As an example, Fig. ^| shows the resulting uncertainty on the quantity 
Y = Xi+X 2 , where the X^ are independent and their uncertainty is described by identical 
asymmetrical triangular distributions. The combined result is obtained analytically using 



1 

1 

\/6 

T3 

72 
1 

L64 
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L96 




0.58. 
0.41. 

: 0.42. 
0.61. 
0.51. 



(21) 
(22) 

(23) 
(24) 
(25) 



E[X] = 0.17 

a(X) = 0.42 

mode = 0.5 

median = 0.23 



E[X] = 0.17 

a{X) = 0.42 

mode = 0.5 

median = 0.23 



E[X] = 0.34 

<r(X) = 0.59 

mode = 0.45 

median = 0.37 





Figure 2: Probability density function resulting from the sum of two quantities, each described 
by an asymmetric triangular p.d.f. with xq = 0.5, Ax+ = 0.5 and Ax_ = 1.5. 



9 ) The so called coherent bet, a recognised normative tool to assess probability, can help a lot to elicitate 
model and parameters of type-B uncertainty. For a concise introduction on its concept see Ref. |To| . An 
extensive discussion about its role to force people to assess uncertainty in non-standardised situations, 
see Ref. [pi. 
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Eq. (|8|). One can see how good the Gaussian approximation already is and how biased a 
result can be, if the best estimate of the sum is performed using mode or median, and if 
the final uncertainty is evaluated with ad hoc rules of the kind shown in the introduction. 

5 Small deviations from linearity 

Let us consider now nonlinearity effects, which are mostly responsible for the published 
asymmetric uncertainties due to systematics. Nonlinearity makes in fact A +ii and A „ 
differ considerably. We treat here only second-order effects. Figure ^| shows an example 
of the transformation of some important p.d.f.'s, all characterized by ELY] = and 
a(X) = 1, while Fig. ^ shows the convolution of the original and of the transformed 
quantities of Fig. |3|. One can see that indeed the p.d.f. of the sum of of both the original 
and the transformed quantities can be described by a Gaussian for the practical purposes 
of interest in uncertainty evaluations. 

In order to simplify the formulae, let us consider first the case of only one input 
quantity and one output quantity (see Appendix B for the general case). Taking the 
second-order expansion, we have 

BY 1 f) 2 Y 

Y = Y(E[X}) + Qx&- E[X]) + 2BX~ 2 {X ~ E[X])2 ■ (26) 

Expected value and variance of Y are then 

E[Y] = Y(E[X}) + ±^f 2 a\X), (27) 

+ \ t E[(x - E[x])4] - a4(x) ] • (28) 

These formulae can be transformed into more practical ones if the derivatives are replaced 
by their numerical evaluations from the ±1 a of X around ELY], which produce variations 
A± in Y. The approximate derivatives evaluated in E\X\ are 

dY 

BX 
d 2 Y 
BX 2 

The formula of the variance 
(£), defined as S = E [(X 
tively. We get finally 

E[Y] 
AY) 

where 5 is the semi- difference of the two shifts [5 = (A + — A_)/2] and A is their av- 
erage [A = (A + + A_)/2]. The interpretation of Eq. (|3lD is simple and corresponds to 
a procedure that some might have already guessed: Asymmetric uncertainties produce a 



2 \a(X) a(X) J 2a(X) ' 
1 ( A+ A_ \ A + - A_ 



a{X) \a{X) a(X)J a 2 (X) 



(30) 



Eq. (^3), can be simplified using skewness (S) and kurtosis 
- E[X]f] /a 3 (X) and K = E [(X - ELY]) 4 ] /a 4 (A), respec- 



Y(E[X}) + 5, (31) 
A 2 + 2A-5 ■ S(x) +5 2 -[IC(X) -1] , (32) 
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Figure 4: Probability density functions of the sum of the quantities Xi and of their nonlinear 
transformations Yi defined in Fig. [3|. 



shift in the best estimate of the quantities. In the case that the dependence between Y 
and X is linear, 5 is ~ and we recover the result given in Section |3[ Note also that the 
second term of Eq. (|32| ) disappears if the distribution describing the uncertainty on X 
is symmetric around E[X], and that the third term plays a second-order role, since the 
difference between A + and A_ is usually smaller than their sum, and JC(X) is around 2 
or 3 for the distributions of interest (see Fig. [I]). 

The extension to several independent input quantities is straightforward, as one only 
needs to add together the individual contributions to expected value and the variance. 
Considering the most common case in which the second and third terms of the r.h.s. of 
Eq. (|3^) are negligible! 10 ^, we obtain the following simple practical formulae: 

E[Y) » F(E[X])+5>, (33) 

i 

a\Y) » £a*. ( 34 ) 

i 

Averaging positive and negative deviations is indeed a good practice, but the shift of 
the central value should not be neglected. For the separation of input quantities into 



fj, r . and influence factor, see Eqs. (p!9|)-(|20|). The formulae for the more general case of 
several output quantities and of correlations among input quantities will be considered in 
Appendix B. 



6 Numerical examples 

Let us go back to the numerical example of the introduction. Those numbers were obtained 
from a quadratic dependence of Y on the influence quantities, each having a slightly 
different functional form and a different model to describe its uncertainty. Including also 
/i r as Xq, we can write the dependence of Y on Xi in the following explicit form: 

3 

Y = Y,a i X l + (3 i Xf } (35) 

where a« and $ are given in Table [I], in which also the uncertainty model is indicated. 
As stated in the introduction, the expression 'reasonable variation of the parameters' was 

10 ) For symmetric distributions the skewness is zero, while the kurtosis is around 3 for the distributions 
of interest and enters with S 2 . 
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Table 1: Parameters of the input quantities used in the numerical example of the text. Xq is 
identified with the value fi r , obtained when X1-3 are equal to their expected values. 



Interpretation 1: 'reasonable variations' = ±1 a for all X; 


Input/Output 


Distribution 


E[X] 


a(X) 


a 


P 


A_ A+ 


X (= Hr) 


Gaussian 


1 


0.05 


1 





+0.050 +0.050 


x x 


Gaussian 





0.3 


0.25 


-0.167 


+0.090 +0.060 


x 2 


Triangular [—1,1] 





0.41 


0.30 


-0.147 


+0.147 +0.098 


x 3 


Uniform [—1,1] 





0.58 


0.225 


-0.078 


+0.156 +0.104 


Y 


~ Gaussian 


0.93 


0.20 








Interpretation 2: 'reasonable variations' 


= ±1 a for fi r 


and Xi] 


±Ax for others 


Input / Output 


Distribution 


E[X] 


a(X) 


a 


(3 


A_ A+ 














(rescaled at 1 a) 


X (= fx r ) 


Gaussian 


1 


0.05 


1 





+0.050 +0.050 




Gaussian 





0.3 


0.25 


-0.167 


+0.090 +0.060 


x 2 


Triangular [—1,1] 





0.41 


0.123 


-0.0245 


+0.054 +0.046 


x 3 


Uniform [—1,1] 





0.58 


0.130 


-0.026 


+0.084 +0.066 


Y 


~ Gaussian 


0.97 


0.13 









intentionally left vague. We consider the two cases in which the variations of non-Gaussian 
quantities correspond to ± 1 a or to ± half-interval, respectively ('interpretation 1' and 
'interpretation 2' in Table [l]). The details of the first evaluation are (see Appendix B for 
the second case and for the values of central moments of higher order) 

E[Y] = 1.00 + ^5, = 1.00 + (-0.015 - 0.026 - 0.0245) = 0.9345, (36) 

a\Y) = ( r r 2 (r) + a s %(F) = (0.05) 2 + (0.1983) 2 =(0.2046) 2 , (37) 

a result which can be summarized as Y = 0.93(0.20), or Y = 0.93 ± 0.20, although the 
latter expression might be misleading, since it is traditionally used as a 68% probability 
interval, which is exactly true only if the p.d.f. is perfectly Gaussianf""^ (see also comments 
in Ref [I9[). The result (|36|)-(|37|) is in excellent agreement with E[Y] = 0.9344 and 



cr(Y) = 0.2046 obtained directly from the p.d.f. of Y estimated by Monte Carlo with 10 
extractions. In contrast, the result obtained combining separately positive and negative 
deviations in quadrature (see introduction) shows a bias which amounts to 35% of a. 



l > As stated above and shown with several figures, in most of the practical cases the Gaussian approx- 
imation is a good one, even if the models describing the uncertainty of the input quantities are not 
Gaussian or if there are some nonlinearity effects. In summarizing the result with a couple of numbers, 
our preference goes to expected value and standard deviation, because these are the parameters which 
matter in further propagation of uncertainty, and this is the most common use of scientific results. 
In this respect we agree with the recommendations of the ISO Guide [M. In the case that the final 
p.d.f. differs considerably from a Gaussian, is indeed a good practice to provide also mode and median 
of the distribution, as well as probability intervals of interests (see e.g. Refs. [^l| and However, 
it is clear that if, for decision problems, one wishes to assess in the most precise way probability 
intervals, the exact form of the p.d.f. is required. For complex problems this evaluation is performed 
by numerical or Monte Carlo techniques. We give in Appendix B approximate formulae of skewness 
and kurtosis, which give an idea of the deviation of the p.d.f. from a Gaussian. Appendix C shows also 
how to get an idea of the shape of the p.d.f., assuming that it is not 'too' different from a Gaussian. 



14 



Y = Y(X) 


f 




f 














\ 

\ 1.3 
\ / 

- 

N n/c 
\ / 

\ of. 

\ / 

N b ' 

\ r ■ 
\ / 

\ /°-* 


\ 

\ 

- 


1 

1 

1 
1 

1 

1 

1 

/ 

/ 

/ 

/ 

/ 

/ 

/ 

\ / 
\ / 


200 
150 
100 
50 
























-3 -2 -1 




12 3 




0.005 


.01 


.015 


. 02 


0.025 


0.03 



Figure 5: Example of non-monotonic relation between input and output quantity. The left plot 
show the parabolic dependence of Y on X (dashed line) and the Gaussian p.d.f. of X (solid 
line). The right plot shows the p.d.f. of Y. 



7 The non-monotonic case 

Sometimes a variation of ± la of an influence parameter might produce values of Y which 
are both above or both below the value obtained with the reference value. Using the 
notation of this paper, A + and A_ have opposite signs in that case. This result indicates 
that the function is not monotonic, and this situation has to be treated with some care. In 
fact, although the formulae derived in this paper do not depend on whether the functions 
are monotonic or not, the transformed distribution can be very different from those of 
Fig. [| and can bring a large non-Gaussian contribution to the overall distribution. As an 
example, let us consider Fig. || which describes an input quantity normally distributed 
around with a = 0.3, a parabolic dependence of Y on X given by Y = 0.167X 2 (i.e. like 
X% of Table [1], but with the a = and j3 reversed in sign, just for graphical convenience). 
The ±lcr variations are A + = +0.015 and A_ = —0.015, but certainly one would not 
quote as the expected value of Y, nor 0.015 its standard deviation. E[X] being at the 
minimum of the distribution, the p.d.f. of X ends sharply at zero, and is very asymmetric. 
In fact it is easy to recognize in f(y) a scale transformation of the y 2 with one degree 
of freedom, namely Y = 0.015 x xi- Expected values and standard deviation are then 
E[Y] = 0.015 and a(Y) = 0.015 x = 0.021. We can compare the result with what we 
get from Eqs. ©-©: 

f A = f E[Y] = + 0.015 = 0.015 

\ 5 = 0.015 ^ \ a 2 {Y) = + + 0.015 2 x 2 = (0.021) 2 W 

The result is exactly the same, as it should be, since in this example the function is 
parabolic and, therefore, there are no approximations in Eqs. (|3~TD-(j3!|). We see that in 
this case only the quadratic terms appear. Similarly, it would be wrong to consider the 
best estimate of Y as equal 0, with an uncertainty equal to the deviation: The result 
would have a standard deviation smaller by v^2, and the best estimate would have a bias 
of —140% of the reported standard deviation. 
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8 Conclusions 



The issue of uncertainty propagation has been reviewed using the probabilistic approach 
nowadays called Bayesian. Following physicists' intuition, this approach allows probabilis- 
tic reasoning to be applied to not exactly-known values of physical quantities. Therefore, 
probabilistic statements then have a clear meaning and can be propagated to the values 
of related quantities using probability calculus. In particular, uncertainties due to sys- 
tematic effects or theoretical inputs can be included in a conceptually easy way. Indeed, 
there is no qualitative distinction between uncertainty due to random and systematic (or 
theoretical) effects, and in this respect we agree fully with the ISO Guide fL9f . 

Approximate formulae have been derived to calculate the first moments of out- 
put p.d.f., taking into account up to second-order effects in the uncertainty propagation. 
Non-linear effects are responsible, together with the less frequent case of asymmetric 
uncertainties on input quantities, of results given with asymmetric uncertainties, if un- 
certainties are evaluated using ad hoc procedures. We have shown that in most cases of 
practical interest, a combinatorial effect similar to the central limit theorem symmetrizes 
automatically the overall uncertainty. It follows that expected value and standard devia- 
tion of the distribution should be the correct way of presenting the result, as also pointed 
out by metrological authorities [ 19| . In the case that the final distribution is considerably 
non-Gaussian, other position and width indicators of the distribution (and whenever pos- 
sible the p.d.f. itself) should also be provided, but expected value and standard deviation 
should always be given, as they are what mostly matters in further propagation of un- 
certainty. We have also shown that asymmetric uncertainty on input quantities and/or 
non-linear effects produces shifts in the expected value of the output quantity with respect 
to the value calculated for the best estimate of the input quantities. // these shifts are not 
applied, the result can be biased. 

As far as modelling the uncertainty due to systematic effects is concerned, we have 
found particularly helpful the concept of type-B uncertainty, according to the ISO termi- 
nology We have shown that the results are rather stable against reasonable variations 
of models and values of the parameters. On the other hand, sticking to the position that 
there is no way of modelling uncertainty due to systematic effects would lead to impossi- 
bility of providing an overall uncertainty, or to providing this uncertainty using arbitrary 
prescriptions which do not have a clear meaning. In other words, we prefer beliefs assessed 
by people we trust, rather than empty prescriptions, in accordance with the scheme: be- 
liefs in, beliefs out; nothing in, nothing out (see Ref. [jT0[] for a discussion about 'belief 
and 'arbitrariness'). This position is based on the assumption that "it is scientific only to 
say what is more likely and what it is less likely" ||25|| . 

We would like to conclude with some remarks about when to stop making systematic 
checks and adding contributions to the overall uncertainty. As eloquently said p6|, "one 



could correlate the result with the phase of the Moon or the position of Jupiter, and find 
most likely no significant effect, with some uncertainty; but certainly we don't want to 
take care of this uncertainty. " Only contributions which are in principle relevant should 
be considered in the uncertainty evaluation. Even if the effect is 'statistically significant', 
one should try to understand if it can physically influence the result, before including it in 
the analysis. There is certainly some subjectivity in deciding what is relevant and what is 
not, but this is consistent with the spirit that "the quality and the utility of the uncertainty 
quoted for the result of a measurement therefore ultimately depend on the understanding, 
critical analysis, and integrity of those who contribute to the assignment of its value" fl9lj . 
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Appendix A — Derivation of Eq. (§) and comparison with Eq. (|3]) 

Equation (§) is one of the most convenient ways to formulate the problem of p.d.f. prop- 
agation, especially for physicists, who are familiar with the Dirac delta. Better known 
formulae which make use of the Jacobian are easily recovered by making use the proper- 
ties of the delta. The same is true for the convolution formula to obtain the p.d.f. of the 
sum of two independent quantities. In this latter case, the derivation is straightforward: 
f(y) = I f( x i) ■ ffa) ■ Kv - x i~ x 2) dxidx 2 = / f(xi) ■ f(y - xi) dx x . Equation (D is 
also important because it justifies the Monte Carlo estimation of the p.d.f. and clarifies 
the role of the various ingredients which enter the game (see comments in footnote 7). 

Although we think that Eq. @ does not need to be proved, as it can be considered 
itself the formulation of the problem, we give here a formal derivation based on the 
properties of the characteristic function. The characteristic function associated with X is 
defined as (see e.g. [E3|): 



= E [e ux ] = J t 



e ltx f(x)dx, (39) 



from which the p.d.f. can be reobtained as 



/(x) = i j e ~ Ux Mt)dt. (40) 
If we have a function Y = g(X), the p.d.f. of Y can be obtained from the following 



characteristic function 23 



fait) = J e it9{x) f(x)dx. (41) 
This property can be extended to a variable Y depending on many variables, i.e. 

M*)= /e ftfl(x) /(x)dx. (42) 

It follows 

f( y ) = i_ J e~ ity (t) Y {t) dt = J e~ ity dt J e it9{x) f(x) dx = J e - i *(^W)/( x ) dxdt. (43) 
Noting that J e %t ^ x ~ yS> dt = 2n5(x — y) we get finally 

f(y)= ff(x)8{y-g(x))dx, (44) 



which is equivalent to Eq. @. 

Once we have proved Eq. (|8|), there is no need to prove also that this is equivalent 
to Eq. (HD, as the latter comes from a general theorem of probability theory. We show 
their equivalence in a simple case, using the example of the offset uncertainty given in 
Section |[ Using the notation of Eq. (|8]), the correspondence of symbols is Y = fi, hi = z 
(with z = 0) and "data" = x. The p.d.f.'s of our input quantities are 

1 _(nr-^) 2 

f(Hr | x) [= f{n | x, z )} = e 2 ^ , (45) 

\Z2ira r 

f(z) = _^ e "^I, (46) 

V2 7T(T 2 
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while the function that relates output to input quantities is n — /J, r + z . Applying Eq. 
and writing explicitly the integration limits, we have 



f{fi\x) 



oo r+oo 



oo 
+00 



(ur—x) 



2 7T O r 



/ 27T(T 2 



_ z 2 

e 5(fi — fi r — z) dfi r dz (47) 



(n-z-x) 
7^2 



2 7T O r 



1 -4j 

e 2<T z dz 



x 2 7T cr. 



(48) 



which is exactly Eq. 



Appendix B — Useful formulae up to second-order approximation 

In the general case of several input quantities and several output quantities the second- 
order expansion is given by 

dY 4 



Y j = y j (E[X]) + ^^(X J -E[X,]) 

i 

1 d 2 Y 



2 ^ dXidXr, 

Lm 



(X l -E[X l ])(X m -E[X m \). 



Expected value, variances and covariances of Yj are then 

d 2 Y, 



Eft] 
a 2 (Yj) 



^i(E[X]) + i £ — Cov(X, X m ) 



2 dXidX n 

Lm 



E 

l,m 



dYj 
dXi 



8Yj_ 

dXrr, 



E 

Lm,n 



+ ^E 



ax ax m ax n 



Cov(Xi,X n 



E[X^X m X n ] 



h,l,m,n 



d 2 Y 



dX h dXj \dX m dX n 



d 2 Y 



Cov(E„E fcJ 



x ^E[X /l X / X m X n ] - Cov(X fe , X) Cov(X m , X„ 
dYj 8Y k 



«9X «9X„ 



Cov(X,X„ 



+ E 



l,m,n 



dY k 8% dYj d 2 Y k 



+ 



dXdX^dXr, dXidX^dX. 



l,m,r 

[x h : 



02 



E 



X^X m X n 



d 2 Y 



dX t dX h 



h,l,m,n 
X ( E I XhXlX m X n 



d 2 Y k 



dX m dX n 



Cov(X fe ,X)Cov(X m ,X„ 



(49) 



(50) 



(51) 



(52) 



where Xj = Xj — [Xj]. The expressions involving expected values of products of more 
than two Xj's can be evaluated by iterating the procedure described here, as they are 
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functions of X,. An interesting case for the applications is, using the same notation as 
at the end of Section |3|, that in which we have N input quantities of the kind and 
all the remaining ones have the meaning of uncorrelated influence quantities. Using the 
derivatives evaluated numerically and the definition of skewness and kurtosis, we can write 
the previous formulae as 



i>N 

+ 2 £ Kit ■ 8fi ■ S{Xi) +J2 5 %- [*W " !] 

i>N i>N 

Cov(/i ri , (JL rk ) + ^ Aji^ki 

i>N 

+ ^2(Kji 8 hi + A u Sji) SiX,) + J2 SjiSki [K(Xi) ~ 1] , (54) 



(53) 



Gav{Yj,Y k 



i>N 



i>N 



where = (A + ji + A_^)/2 and Sji = (A +J j — A_jf)/2 and A±ji are defined according to 
the convention of footnote |8|. As explained in the text, the terms depending on skewness 
and kurtosis are negligible in most practical cases. 

It might be useful to have also the approximate expressions of skewness and kurtosis 
of the output quantities Yj, to have a rough idea of how much their p.d.f.'s differ from 
the Gaussian. Since the formulae become quite awful, we only consider the most common 
case of uncorrelated input quantities. 



S(Yj. 



3 k • -!]+£<%■ [£{Xi) ~ 3 K(Xi) + 2] 



(55) 



+ 4 w + 6 E 5 l 4 • pw - 2 5 ( x ^)] 

£ k{x,) A- t +6j2 K Kn + 6 E ^ • [*w - ^ 
+ 6 E 4 4 ■ pot - 2 *w + !] + 6 E $ • [*w - 1] • [*w - 1] 

i Z>m 

+ E^- [°w + 6 *w - 4 f w - 3 i + 4 E ^ 4 • [ p ( Xj ) - 5 

i i 

+ 4 E^4 • MXi) -3V(Xi) + 3 SiX,)} 



+ 24 E ^ 7 ^ m ^ l ^OrnSiXi) SiX r , 



l>m 



(56) 



where 



V{X) = E[(A-E[A]) 5 ]/a 5 , £{X) = E[(X - E[X]) e ] /a e , 

H{X) = E[(A-E[A]) 7 ]/^ 7 , 0{X) = E[(X-E[X]) 8 ]/a 8 

are higher order scaled central moments, analogues of skewness and kurtosis. Table |] gives 
these higher order moments for most relevant p.d.f.'s used to assess type-B uncertainty. 
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Table 2: The rescaled central moments of fifth to eighth order for the distribution of Fig. |j. The 



par ameters a, (3 and 7 are given in Eqs. (|57|)-(|59|). 



Distribution 


V(X) 


£{X) 


H{X) 


O(X) 


Gaussian 





15 





105 


Uniform 





27/7 





9 


Triangular 





54/7 





144/5 


Asymmetric triangular 


P/a 5 


2/7 x (31 - 27a) 


7/a 7 


16/5 x (13 -27a) 



The moments of the asymmetric triangular distribution are given in terms of the following 
parameters: 

a = Ax 2 _Ax 2 + (Ax_ + A + ) 2 / (Ax 2 _ + Ax^A + + Ax 2 + ) 3 , (57) 
P = -(Ax b _ + 5/2 AxiAx+ + AxlAx 2 + - (58) 

Ax 2 _Ax 3 + - 5/2 A_A* - A^)/425.2 , 
7 = -(2Aa;L + 7Ax 6 Ax + + 9AxiA^ + 5AxiAa; 3 - (59) 

5 Ax 3 _ Ax\ - 9 Ax 2 _Ax\ - 7 A_ A^ - 2 A 7 _) /2915.0 . 

As a numerical example, continuing on from Section [| we compare the results obtained 
with the approximate formulae with those obtained calculating the moments directly 
from the p.d.f. estimated by Monte Carlo. The overall results are given in Table [3[ The 
agreement is well above that needed for practical purposes. 



Table 3: Comparison between the moments evaluated using approximate formulae and those 
obtained by Monte Carlo, based on the examples of Table |l[ 





'Interpretation 1' 


'Interpretation 2' 




Monte Carlo 


Approx. formulae 


Monte Carlo 


Approx. formulae 


E[Y] 


0.9344 


0.9345 


0.9722 


0.9720 


a(Y) 


0.2046 


0.2046 


0.1297 


0.1295 


son 


-0.370 


-0.372 


-0.318 


-0.321 


£(Y) 


2.857 


2.859 


3.076 


3.082 



Appendix C — Approximate evaluation of the p.d.f. from the first four 
central moments 

We have seen that the propagation of uncertainty can be solved either exactly, at the cost 
of having to compute complicated multidimensional integrals, or by using approximate 
formulae which only give the first moments of the distribution. The probabilistic interpre- 
tation of the latter case is based on the assumption that the final p.d.f. is approximately 
Gaussian, as is often so in most cases of practical interest. We would like to show here an 
intermediate solution to the problem, which allows the shape of the p.d.f. of each Yj to 
be estimated starting from the approximate evaluation of the first four central moments 
(see Appendix B), and solving only one integral in only one variable. 
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0.25 0.5 0.75 1 1.25 1.5 1.75 2 0.25 0.5 0.75 1 1.25 1.5 1.75 2 

Figure 6: The distribution given by the complete uncertainty propagation (dots) is compared 
with the approximate probability density function for interpretations 1 (left side) and 2 (right 
side) . 

Take the generic output variable Y. Its characteristic function </>y(t) can be ex- 
pressed as 

M*) = «pff^(»t) B J , (60) 

where Xn is the nth 'semi-invariant' Semi-invariants can be expressed in terms of 
central moments, the first four being 

Xi = E(X), 

X2 = O* , 

Xs = <x 3 S(X), 

X4 = a 4 [/C(X)-3]. 

An approximate evaluation of the p.d.f. can be obtained under the assumption that the 
first four moments, for which we have obtained approximate expressions, are the most 
relevant ones to characterize the p.d.f. . This implies that the p.d.f. 's of interest do not 
differ 'too much' from a Gaussian. In practice this means that they are unimodal, have 
no discontinuity, and that the probability mass is concentrated around a few standard 
deviations from the expected value. Assuming, then, that all semi-invariants Xn with 
n > 4 are the same as for a Gaussian distribution, namely zero, the characteristic function 
becomes 

Mt) = exp (i t E[Y) - l -t 2 a 2 (Y) - l -t 3 a 3 S(Y) + 1 1 4 a 4 \K{Y) - 3]) . (61) 

Inserting Eq. ( |B~TD in the anti-transformation formula [see Eq. (|4*3|)1 and taking only the 
real part of the solution, since the p.d.f. is real, we get 

f(y) = i- J exp (~aH* + 1 tV 4 [K{Y) - 3]) cos (t [E[Y] - y] - if 8 a 3 5(F)) dt . 

(62) 

Applied to our examples of Section |] and using the values of Table [3], we get the p.d.f. 's 
drawn with the continuous curves of Fig. |6] [the integral ( |62D has been solved with Math- 
ematical. For comparison the dotted curves show the estimations of p.d.f.'s obtained by 
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Monte Carlo. The agreement as it can been judged by eye is excellent. However, since 
the method is approximate, there are slight problems of normalization and positiveness. 
But these problems affect the tails of the distribution and are not really relevant if one 
is interested in having an idea of the shape under the assumptions of the approximation. 
Note also the divergent term in Eq. (|62|) for /C > 3. But in practical cases the kurtosis 
is never much larger than this value. In fact one starts usually from distributions which 
have JC(X) < 3 (see Fig. HD and, thanks to the central limit theorem, there is a natural 
tendency to have JC(Y) ~ 3. Therefore, in case of values of kurtosis slightly larger that 
3, a good approximation is to limit it at 3. This approximation has been indeed applied 
to the 'interpretation 2' of example of Table |3], and the resulting p.d.f. is still in excellent 
agreement with the Monte Carlo evaluation (see right hand plot of Fig. [|). 
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